/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org.

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Functionality for reading TIPSY files. \ingroup sfrhist

// $Id$

#include <iostream>
#include <string>
#include <iterator>
#include <algorithm>
#include "preferences.h"
#include "snapshot.h"
#include "constants.h"
#include "boost/lambda/lambda.hpp"
#include "misc.h"

using namespace std;

template <typename input_iterator, typename output_iterator>
output_iterator copy_n (input_iterator& source, size_t n,
                         output_iterator destination)
{
  while (n--) {
    *destination++ = *source++;
  }
  return destination;
}

template <typename input_stream, typename output_iterator>
output_iterator read_n (input_stream& source, size_t n,
                         output_iterator destination)
{
  while (n--) {
    double v;
    source >> v;
    *destination++ = v;
  }
  return destination;
}

template <typename input_stream>
void skip(input_stream& i, int n)
{
  double d;
  while(n--) i >> d;
}

/** Loads a TIPSY snapshot file by determining if it's binary or ascii
    and calling the appropriate function.  */
bool mcrx::Snapshot::loadtipsy(const string& snapfn)
/** Loads a Tipsy snapshot file */
{
  ifstream snap;
  snap.open(snapfn.c_str());
  if(!snap) {
    cerr << " Error opening snapshot file: " << snapfn << endl;
    return false;
  }
  int junk;
  snap >> junk;
  snap.close();
  if((junk>0) &&(junk<1000000000)) //  what is the heuristic here?
    return loadtipsy_ascii(snapfn);
  else
    return loadtipsy_bin(snapfn);
}

/** Loads a Tipsy ASCII snapshot file */
bool mcrx::Snapshot::loadtipsy_ascii(const string& snapfn)
{
  using namespace boost::lambda;

  name_ = snapfn;

  ifstream snap;
  snap.open(snapfn.c_str()); //Check that this will work
  if(!snap) { 
    cerr << " Error opening SNAP file!\n";
    return false;
  }

  cerr << "Reading TIPSY ASCII snapshot file " << snapfn << endl;

  // numbers of particles -- No idea if this is needed.  We'll see
  int n_species = 5;
  int n_tot;
  int nsph;
  int nstar;
  const int ndisk = 0;
  const int nbulge = 0;
  snap >> n_tot;        //changed this from _ntot to n_tot, hope that is ok
  snap >> nsph >> nstar;
  const int nhalo = n_tot - nsph - nstar ;

  // dimensions is always three, and time is actually expansion
  // factor, so we demand the time be set in preferences
  // (in internal units)
  double junk;
  snap >> junk >>  a_;
  time_ = prefs_.getValue ("snaptime",double());  

  const bool comoving = prefs_.defined("comoving_units") ?
    prefs_.getValue("comoving_units",bool()) : false;
  if (!comoving) {
    // if not comoving, we simply ignore the a value.
    a_=1;
    cout << "Comoving units are not used" << endl;
  }
  else
    cout<<"Comoving units are used, expansion factor is a="<< a_ << endl;

  // No need to check for BHs

  // length in kpc
  lengthunit_ = a_*(prefs_.getValue("UnitLength_in_cm",double()))/    
    (constants::kpc*1e2);
  // mass in Msun
  massunit_ = (prefs_.getValue("UnitMass_in_g",double()))/
    (constants::Msun*1e3);
  // time in yr
  timeunit_ = prefs_.getValue("UnitLength_in_cm",double())/(prefs_.getValue("UnitVelocity_in_cm_per_s",double())*constants::yr);  

  G_=0; // Note: _G will be incorrect until makephysicalunits() is run, we can figure out a _G here and we should  

  cout<<"Total Number of particles: "<<n_tot<<endl;
  cout<<"Number of Gas: "<<nsph<<endl;
  cout<<"Number of stars: "<<nstar<<endl;
  cout<<"Number of halo: "<<nhalo<<endl;

  // because the tipsy file does not have ID numbers, we don't have to
  // deal with any reordering so we can read directly into the
  // vectors.
  
  // Read masses
  read_n(snap,nsph,back_inserter(gas_particles.m_)); //gas 
  read_n(snap,nhalo,back_inserter(halo_particles.m_)); //halo
  read_n(snap,nstar,back_inserter(nstar_particles.m_)); //star
  assert(gas_particles.m_.size()==nsph);
  assert(halo_particles.m_.size()==nhalo);
  assert(nstar_particles.m_.size()==nstar);

  // Read x..y..z.. and vx vy vz.  since they have to be assembled
  // into three-vectors, we read them into temporary vectors first
  // fields are ordered halo, gas, nstars
  vector<double> tx,ty,tz,tvx,tvy,tvz;
  read_n(snap,n_tot,back_inserter(tx));
  read_n(snap,n_tot,back_inserter(ty));
  read_n(snap,n_tot,back_inserter(tz));
  read_n(snap,n_tot,back_inserter(tvx));
  read_n(snap,n_tot,back_inserter(tvy));
  read_n(snap,n_tot,back_inserter(tvz));

  // now assemble
  // gas particles are first
  gas_particles.pos_.resize(nsph);
  gas_particles.vel_.resize(nsph);
  for(int i=0;i<nsph;i++) {
    gas_particles.pos_[i] = vec3d(tx[i],ty[i],tz[i]);
    gas_particles.vel_[i] = vec3d(tvx[i],tvy[i],tvz[i]);
  }

  // then halo
  halo_particles.pos_.resize(nhalo);
  halo_particles.vel_.resize(nhalo);
  for(int i=0;i<nhalo;i++) {
    halo_particles.pos_[i] = vec3d(tx[i+nsph],ty[i+nsph],tz[i+nsph]);
    halo_particles.vel_[i] = vec3d(tvx[i+nsph],tvy[i+nsph],tvz[i+nsph]);
  }
  // finally stars
  nstar_particles.pos_.resize(nstar);
  nstar_particles.vel_.resize(nstar);
  for(int i=0;i<nstar;i++) {
    nstar_particles.pos_[i] = vec3d(tx[i+nsph+nhalo],ty[i+nsph+nhalo],
				    tz[i+nsph+nhalo]);
    nstar_particles.vel_[i] = vec3d(tvx[i+nsph+nhalo],tvy[i+nsph+nhalo],
				    tvz[i+nsph+nhalo]);
  }

  //  cout<<"Postions and Velocities Read"<<endl;

  // halo softening length
  // This is never comoving.
  read_n(snap,nhalo,back_inserter(halo_particles.h_));
  transform(halo_particles.h_.begin(), halo_particles.h_.end(),
	    halo_particles.h_.begin(),
	    _1 * (1.0/a_) );

  // stellar softening length goes into _h for the new stars
  read_n(snap,nstar,back_inserter(nstar_particles.h_));

  // We need to back-convert the expansion factor, because the
  // softening lengths are physical in the file, not comoving.
  // (star_radius_factor is applied in the sfrhist_snapshot
  // constructor.)
  transform(nstar_particles.h_.begin(), nstar_particles.h_.end(),
	    nstar_particles.h_.begin(),
	    _1 * (1.0/a_));

  //Will wait a while to do the the gas softening as we need to read
  //it in from a different file

  vector<double> trho;  //Not sure that density is written, keep it for now
  read_n(snap,nsph,back_inserter(trho));

  vector<double> temps;
  read_n(snap,nsph,back_inserter(temps));

  // Gasoline feedback does not make a distinction between thermal and
  // effective temperature, so they are the same.
  gas_particles.temp_.assign(temps.begin(),temps.end());
  gas_particles.teff_.assign(temps.begin(),temps.end());

  cout<<"Temperature read"<<endl;

  // now smoothing lengths and copy temporary into real vector
  // actually they say this is the gravitational softening length, so
  // instead we read the smoothing lengths from the supplied file
  // the first nsph numbers should be the gas particles
  if(!prefs_.defined("smoothing_file" )) {
    cerr << "Error: no smoothing lengths file specified.  Add keyword \"smoothing_file\"" << endl;
    return false;
  }
  const string smooth_file_name =  
    word_expand(prefs_.getValue("smoothing_file",string()))[0];
  cout << "Reading smoothing lengths from file " 
       << smooth_file_name << endl;
  ifstream smoothfile(smooth_file_name.c_str());
  if(!smoothfile) {
    cerr << "Could not open smoothing lengths file " << smooth_file_name << endl;
    return false;
  }
  int n_smooth;
  smoothfile >> n_smooth;
  assert(n_smooth==n_tot);
  read_n(smoothfile, nsph, back_inserter(gas_particles.h_));

  // and remove softenings from the main file
  skip(snap,nsph);

  // metallicities (first gas then stars)
  read_n(snap,nsph,back_inserter(gas_particles.z_));
  read_n(snap,nstar,back_inserter(nstar_particles.z_));
 
  // formation time
  read_n(snap,nstar,back_inserter(nstar_particles.tf_));
 
  //potential is unneeded
  skip(snap,nsph+nhalo+nstar);

  gas_particles.id_.resize(nsph);
  halo_particles.id_.resize(nhalo);
  nstar_particles.id_.resize(nstar);
  nstar_particles.pid_.resize(nstar);
  // make up ID numbers and order maps out of thin air.
  // *** THIS MEANS WE CAN NOT CORRELATE IDS BETWEEN SNAPSHOTS ***
  for (int i=0; i<nsph;++i) {
    gas_particles.id_[i] = i;
    gas_particles.order_[i]=i;
  }
   for (int i=0; i<nhalo;++i) {
    halo_particles.id_[i] = i+nsph;
    halo_particles.order_[i+nsph]=i;
  }
   for (int i=0; i<nstar;++i) {
     // we have no parent ID info so we point it to self (this will
     // cause failures if you expect parent particles to be gas)
    nstar_particles.id_[i] = i+nsph+nhalo;
    nstar_particles.pid_[i] = i+nsph+nhalo;
    nstar_particles.order_[i+nsph+nhalo]=i;
  }
   // cout<<"Order Set"<<endl;
 
  // sfr doesn't exist
  gas_particles.sfr_.assign(nsph,0.);

  // creation mass of Star particles is from preferences
  // (gas particles have just current mass)
  nstar_particles.mcr_.assign(nstar,prefs_.getValue("m_star_creation",double()));

  snap.close();

  cout<<"File closed"<<endl;

  // *** IF THERE ARE BLACK HOLE PARTICLES THEY NEED TO BE MOVED TO BH PARTICLE SET ***

  return true;
}


/** Loads a binary TIPSY snapshot file. This function has not yet been updated to the new format.  */
bool mcrx::Snapshot::loadtipsy_bin (const string& snapfn)
{
  cerr << "Reading binary tipsy files is not implemented.\n";
  return false;
}

